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ABSTRACT 

We present the results from three dimensional hydrodynamical simulations of decay- 
ing high-speed turbulence in dense molecular clouds. We compare our results, which 
include a detailed cooling function, molecular hydrogen chemistry and a limited C and 
O chemistry, to those previously obtained for decaying isothermal turbulence. 

After an initial phase of shock formation, power-law decay regimes are uncov- 
ered, as in the isothermal case. We find that the turbulence decays faster than in the 
isothermal case because the average Mach number remains higher, due to the radiative 
cooling. The total thermal energy, initially raised by the introduction of turbulence, 
decays only a little slower than the kinetic energy. 

We discover that molecule reformation, as the fast turbulence decays, is several 
times faster than that predicted for a non-turbulent medium. This is caused by moder- 
ate speed shocks which sweep through a large fraction of the volume, compressing the 
gas and dust. Through reformation, the molecular density and molecular column ap- 
pear as complex patterns of filaments, clumps and some diffuse structure. In contrast, 
the molecular fraction has a wider distribution of highly distorted clumps and copious 
diffuse structure, so that density and molecular density are almost identically dis- 
tributed during the reformation phase. We conclude that molecules form in swept-up 
clumps but effectively mix throughout via subsequent expansions and compressions. 

Key words: hydrodynamics - shock waves - turbulence - molecular processes - ISM 
- clouds - ISM: kinematics and dynamics 



1 INTRODUCTION 



tance t< 



An unc uisUiiding uf ImbuluiiLU it uf fmidamuiiUl nnpnr 



models for molecular clouds (e. g. Vazqucz-Scmadcni et al. 
(199rj); Mac Low et al. (199&t ); |Padoan et al. (1998j ); |5toiK 



many research areas in astrophysics, hi particular, 

compressible turbulence plays a central role in the process 
by which st a rs are formed (review ed by Vazquez- Semadeni 
et al. (200C| ); Padoan el al. (2001 )). Stars form ollL uf dimm 1 
clouds bf molecular gag, which have formed out of diffuse 
clouds of atomic gas. The turbulent energy deduced from 
observations of the molecular gas is sufficient to delay the 
gravitational collapse, making molecular turbulence of great 
interest. The purpose of this study is to provide insight 
through the first three-dimensional simulations of molecu- 
lar turbulence. 

Isothermal or polytropic equations of state have been 
the central premise of previous investigations of turbulent 
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et al. (199S )). This simplification has the advantage that pa- 



rameter space can be explored with numerical simulations 
and the influence of magnetic fields and gravity ca n be de- 



termined in s ome detail (e.g . Balsara et al. (2001 ) 
et al. (1995[ ); |Klessen (200(| ; [Hcitsch et al. (200]] ) 



Passot 



The 

isothermal approximation is indeed valid in many specific 
models where molecular cooling is efficient and low gas tem- 
peratures are maintained. Provided the temperature is low, 
one can argue that an isothermal equation of state is as 
good as any other, with little feedback from the the thermal 
pressure on the dynamics. 

However, for cloud formation and evolution molecular 



Lim 



chem i stry and cooling is critical ( Langer et al. 200C ; 
2001 ; Lim et al. 1999 ). Molecular hydrogen forms most ef- 



ficiently where the gas is warm but the grains are cool (Eb 
forms mainly when atoms combine after colliding and stick- 
ing to dust grains). Simple molecules like OH, CO and H2O 



2 Georgi Pavlovski et al. 



form in the gas phase with H2 as the reactive agent. These 
molecules are not only important coolants, but associated 
emission lines provide a means of measuring the cloud prop- 
erties. Molecules are dissociated as a conse quence of fas t 
shocks, UV radiation, X-rays and cosmic rays (Herbst 2000). 
We thus need to study molecular turbulence to determine 
the distribution and abundances of molecular species. 

Even in cool optically thin regions, as assumed here, 
molecular pressures may influence the dynamics. Strong 
shock waves may be driven from hot atomic regions, formed 
by shock waves, into the proximity of cold molecular gas. 
Areas devoid of dust, or in which molecule formation is in- 
efficient (e.g. because atoms will evaporate off warm grains 
rather than recombine) may attain and maintain pressures 
as high as turbulent pressures. 

The rate of decay of supersonic turbulence is important 
to the theory of molecular clouds. A possible consequence 
of the rapid decay of kinetic energy is that the turbulent 
clouds we observe have short lives. The short lives, however, 
may not provide sufficient time for relaxation into thermo- 
dynamic and chemical equilibrium. Hence, cloud structure 
and content evolve simultaneously as a cloud evolves dy- 
namically, rather than proceeding through a series of equi- 
librium states. Therefore, conclusions based on assuming a 
molecular fraction to be a func tion of density alone (e.g. 
Ballesteros-Paredes et al. (1999)) might not always be ac- 
curate. We note, however, that the conclusion reached by 
Ballesteros-Paredes et al. (1999) of rapid cloud formation is 
supported by this study. 

Decaying supersonic turbulence is the subject of this 
initial study. The work is a direct extension, in terms of code 
an d initial conditions, o f isothermal simulations presented 
by Mac Low et al. (199£). Our goals here are as follows: 



• To directly compare the decay of th e kinetic energy 
Eki n wi th those derived in the papers of (p ac Low et al 



199 q) a nd ( Smith, Mac Low fc Zuev 2000| ; pmith, Mac Low 



Heit ach 200C| ). 

• To determine the distribution of the molecular column 



density as opposed to th e total column density (Vazquez- 
Semadcni & Garcia 2001) 



• To study the rate of formation of molecules. 

We omit in this study magnetic field, self-gravity and pho- 
todissociating radiation. We begin with a fully molecular 
cube of dense gas and apply an initial turbulent field of ve- 
locity perturbations. The turbulence was chosen to be suf- 
ficient to dissociate the gas in one case (root mean square 
(rms) speed of GOkms" 1 ) and to leave the gas molecular in 
another case (rms speed of lSkms^ 1 ). 



2 METHODS 



2.1 The equations 

Numerically we solve the time-dependent flow equations: 

dp 



ot 



^>+<».v>„ = -iv P , 



<2i 
dt 



+ v ■ V e = -pV • v + A (T, n, /) , 



(1) 
(2) 
(3) 



djfn) 
dt 



+ V-{fnv)=R(T,n,f)-D(T,n,f), (4) 



where n is the hydrogen nuclei density, e is the internal en- 
ergy density and / is the molecular hydrogen abundance (i.e. 
71h 2 = fn). We consider the g £LS £LS Et mixture of atomic and 
molecular hydrogen with 10% of helium (i.e. hhc = 0.1 n). 
Therefore, the total particle density is n to t = (1.1 — /) re, 
mass density is p = lAnmu and the temperature is T = 
p I (k ntot). An internal energy loss term has been added to 
the r.h.s. of the energy flow equation: A is the loss of en- 
ergy through radiation and chemistry per unit volume. The 
function consists of 11 distinct parts, summarised below. R 
and D are reformation and dissociation rates of molecular 
hydrogen respectively. A detailed descriptio n of the chem- 
istry can be found in Smith fc Rosen (2002). Formulae for 



en (20( 
dixg 



A, R and D are given in the Appen 



2.2 The numerical model: ZEUS-3D 



1992a 



code. 



As a b asis, we employ the ZEUS-3D code ( Stone fc Norman 
b|). This is a second-order, Eulerian, 



finite-difference 

We study here compressible hydrodynamics without 
gravity, self-gravity or thermal conduction. No physical vis- 
cosity is modelled, but numerical viscosity remains present, 
and a von Neumann artificial viscosity determines the dissi- 
pation in the shock front. 

Further coding for molecular chemistry and molecular 
and atomic cooling has been added. Our ultimate goal is 
to develop a reliable code with which we can tackle three 
dimensional molecular dynamics, later adding self- gravity, 
magnetic fields, ambipolar diffusion and radiation. We thus 
restrict the cooling and chemistry lists to just those items 
essential to the dynamics. We have employed the simu lta- 
neous implicit method discussed by Suttner et al. (1997) in 
which the time step is adjusted so as to limit the change in 
internal energy in any zone to 30%. 

The cooling is appropriate for dense cloud material 
of any atomic- molecular hydrogen mixture. We include H2 
ro-vibrational and dissociative cooling, CO and H2O ro- 
vibrational cooling, gas-grain, thermal bremsstrahlu ng and 



a steady-state approximation to atomic cooling (see Smith 
h Rosen 2002 , Appendix A) and Appendix |X| of the article 



We take a very basic network of chemical reactions. 
Time-dependent hydrogen chemistry is included, but C and 
O chemistry is limited to the reac tions with H and H2 w hich 
generate OH, CO and H 2 (see |5mith fc Rosen 2002| , Ap- 
pendix B) and Appendix A| of the article. Equilibrium abun- 
dances are calculated, which are generally accurate within 
the shocks where molecules are rapidly formed and de- 
stroyed. In cold molecular gas, however, the available oxygen 
will probably remain primarily in the form of water even if 
it has not been shock-processed. Therefore, the code is con- 
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Table 1. Power law indices r) of kinetic energy decay. The rj 
coefficients are derived from linear fitting on the time interval 
[90,600] years. Average error ±0.006. 



resolution 



initial rms velocity for the run 



15 km s 30 km s 60 km s 
-1.33 -1.33 
-1.38 -1.36 
-1.48 -1.52 
-1.51 -1.60 

interval [90,300] years; the run 



32 3 -1.19 
64 3 -1.32 
128 3 -1.45 
256 3 -1.34* 

- linear fitting on time 



terminated half way due to numerical difficulties 



structed to follow the shock-enhanced chemistry and cooling 
rather than the cold molecular gas. This will generally not 
be a problem at the high densities where CO cooling and 
gas-grain heating and cooling determine the properties of 
the low-temperature gas. 

The formula for H2 formation is critical to the results. 
Reformation takes place on grain sur faces with the rate 
taken from (iHollenbach & McKee 19781): 



k R = (3 x KT 18 [cmV 11 



r o. 5/a 



1 + 0.04 (T + T g ) - 5 + 2 x IO- 3 T + 8 x 10- e T 2 
with a factor 



fa 



10000 exp(-600/T s ))" 



(•») 



(6) 



which means, with T g — 20 K assumed, that f a is quite 
close to unity in our simulations. Recen t experimental re- 



immariood by Katss ot al, (1900), indicate that the 



suits, 

H_> formation rate could be considerably lower — The inher- 



ent grain mass and uniform space distribution adopted here 
also influence the results. This would imply that the abso- 
lute reformation times derived here may need revision; we 
intend to explore a range of possibilities in a following study. 

The computations were performed on a Cartesian grid 
with uniform initial density and periodic boundary condi- 
tions in every direction, to simulate a region internal to a 
molecular cloud. The hydrogen was fixed to be fully molecu- 
lar at the beginning (fraction: / = 0.5). Initial stress was in- 
troduced by Gaussian perturbations applied to model veloc- 
ities with a flat spectrum extending over a narrow range of 
wave numbers given by 3 < |fe| < 4. We have performed runs 
with root mean square velocities of 15 km s _1 , 30 km s _1 
and 60 km s~ 
128 3 and 256 



and computational domain sizes of 32 3 , 64 3 , 
zones. The number density has been taken 



to be n = 10 cm , physical box size L — 10 



16 



and 



the initial temperature has a homogeneous distribution of 
Ti — 100 K. Abundances for helium, free oxygen and car- 
bon were taken as 0.1, 5 x 10 -4 and 2 x 10 -4 . 



3 DATA ANALYSIS 
3.1 Decay rates 

The decay in total kinetic energy is displayed in Fig. [j]. An 
initial phase of evolution during which waves steepen and 
velocity gradients transform into shocks is apparent. The 



Table 2. Power law indices fi of thermal energy decay. The fi 
coefficients are derived from linear fitting performed on the time 
interval [90, 600] years. Average error ±0.008. 



resolution 



initial rms velocity for the run 



32 3 
64 3 
128 3 
256 3 



15 km s 1 30 km s 1 60 km s 1 
-0.78 -0.82 -0.47 1 " 

-0.97 -0.99 -0.34* 

-0.82 -0.81 -0.63+ 

-0.85* -0.75 -0.62t 

* - linear fitting on time interval [90, 300] years; the run termi- 
nated half way due to numerical difficulties 

T - linear fitting on the time interval [140, 600] years, present final 
decay behaviour after initial strong cooling regime (see Fig. ^) 

Table 3. Power law indices v of Mach number decay. The v 
coefficients are derived from linear fitting performed on the time 
interval [90, 600] years. Average error ±0.005. 



resolution 



initial rms velocity for the run 



15 km s 1 30 km s 1 60 km s 1 
-0.08 -0.13 -0.46 1 " 

-0.16 -0.15 -0.55* 

-0.21 -0.16 -0.49t 

-0.23* -0.24 -0.54t 

* - linear fitting on time interval [90, 300] years; the run termi- 
nated half way due to numerical difficulties 

' - linear fitting on the time interval [140, 600] years, present final 
decay behaviour after initial strong cooling regime (see Fig 



32 3 
64 3 
128 3 
256 3 



shock formation time is short, with a timescale 

L ._ „ L 10[kms _1 ] 3.5 



ti 



2 km V 



= 45.3 x 



in u rms 



10 16 [cm] 



k',, 



M. (7) 



where fc m is de f ined as the initial mean wavenumb er (Mac 



Low et al. 1998 



Stone et al. 1998| ; |Smith et al. 2000[ ). Hence, 
shocks develop faster in the 60 km s _1 turbulent field, as 
evidenced by the immediate power-law decay. 



if), average molecular fraction 
: 1 5 km s" 1 

: 30 km s"' 




300 
years 



Figure 4. Average molecular fraction as a function of time. There 
is virtually no dissociation during the run with an initial rms of 
15 km s _1 - all four resolution curves merge in one straight line 
to the 0.5 level. 30 km s _1 - intermediate case; 60 km s _1 - 
practically all molecules are temporarily destroyed. 
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E kin decay, rms velocity 15 kms 



E kin decay, rms velocity 30 kms 




1.5 2.0 2.5 

og(time[yr]) 

E kin decay, rms velocity 60 kms" 




1.5 2.0 2.5 

bg(time[yr]) 




1.5 2.0 2.5 

bg(time[yr]) 



resolution 32 3 

resolution 64 3 

resolution 1 28 3 

resolution 256 3 



Figure 1. Average kinetic energy 
as a function of time; log-log plot. 



E m decay, rms velocity 15 kms 



E lh decay, rms velocity 30 kms 



5 7 

r 
5 : 

o : 
s : 
o ; 



1.5 2.0 2.5 

og(time[yr]) 



El,, decay, rms velocity 60 kms 




1.5 2.0 2.5 

og(time[yr]) 




2.0 2.5 
og(time[yr]) 



resolution 32 

resolution 64 3 

resolution 1 28 3 

resolution 256 3 



Figure 2. Average thermal en- 
ergy as a function of time; log-log 
plot. 
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A rapid power-law decay of kin etic energy with time 



was un covered in the isothermal cas e (3tone et al. 1998 
1998| 



Mac 



Low et |al. 1998| ; |Smith et al. 20001 ). This actually followed 
an initial phase of evolution during which velocity gradients 
increase and shocks appear. The time scale for decay is of the 
order of a wave crossing time, which implied that turbulent 
molecular clouds must be short lived if the turbulence is not 
regenerated. 

The first result here is the confirmation that molecu- 
lar turbulence behaves similar to isothermal turbulence, but 
with higher rates of decay. The curves of average kinetic en- 
ergy decay are shown in Fig. |l| Power laws deduced from 
linear fitting (log Ek in oc r\ log t) of the curves are shown in 
Table 0. 

We have found that the formula 



Ekin = -£"0(1 + — ) V , 
ti 



(8) 



fits the data quite well. However, values of \r)\ derived from 
from fitting the curve ([]) are slightly larger due to nonzero 
values of ti . This functional dependence explains why values 
of the average kinetic energy at the end are compatible, 
while at the beginning they are different by the factor of 4. 

The power-law index r\ increases with v rms . This is con- 
sistent with the high power-law index fo und for high Mac h 
number, M — 50, isothermal turbulence (Smith et al. 2000). 

The second result is that the total thermal energy 
(Eth = nfcT(3.3 — /)/2) also decays with time , as opposed to 
the imposed constant value for the isothermal case. The av- 
erage thermal energy density behaviour is presented in Fig. ^ 
and the corresponding decay rates fi: Eth oc t M are shown 
in Table |j| The thermal energy aids our description of the 
physics. We have found that it also has a power-law decay 
for the low speed case in which molecules are not destroyed. 
The thermal energy decays, although not as fast as the ki- 
netic energy. Thus the Mach number remains high relative to 
the isothermal case, which may well account for the steeper 
kinetic energy decay rate. Eventually, however, the gas cools 
down to just below the grain temperature and the turbulent 
heating input falls off. At this temperature range, however, 
the numerical code is not accurate enough to resolve the 
minimum. 

The fact that the thermal energy exhibits a dependence 
on numerical resolution (to a higher degree than kinetic en- 
ergy does) is due to the importance of cooling layer res- 
olution for the energy balance and smoothening effects on 
coarser grids. This implies lower values of thermal energy on 
finer grids, although a smoothening effect for the 64 3 run led 
to higher (than the 32 3 run) values of the average thermal 
energy. 

With an initial rms speed of 60 kms -1 , the molecules 
predominantly dissociate within 40 years. Cooling associ- 
ated with the dissociation and shock heating are competing 
processes (Fig. ^). After 40 years, there is a steep decay in 
thermal energy as trace CO and H2O molecules efficiently 
cool gas between 100 K and 8000 K. As confirmation, we 
note that the number of zones within which the tempera- 
ture is very high (> 10 3 K) falls from a percentage « 45% 
to ~ 1% during the first 150 yr. The gas which had become 
almost fully atomic (see plots of average molecular fraction 
in Fig. ^ and histograms on Fig. [5]) , is again over 60% molec- 
ular hydrogen after 300 yr. As the molecules reform, the gas 



is heated since the molecule energy released on reformation 
is channelled collisionally into the gas (rather than being ra- 
diated) at the high density of these simulations. Therefore, 
the thermal energy subsequently decays slower than in all 
other cases, approaching a shallow power-law value. These 
changes are not directly reflected in the kinetic energy; how- 
ever, the rapid gas cooling maintains a high Mach number, 
which leads to the fast decay of kinetic energy. 

For reference, the energies are plotted in erg cm -3 , with 
the initial kinetic energy given by 



Ei = \pv 2 rms = 1.2 x 10 6 . 

2 \ 10 b [cm ^ 



10 [kr 



[ergc 



(9) 



and the final (thermal) energy, after complete molecule ref- 
ormation (7 = 10/7), is 



E f 



1 -nkT= 1.6 x 10" 9 n ° 



1-7 



10 6 [cm" 3 
Tf 



10 [K] 



erg cm 



(10) 



Plots of the average Mach number are shown in Fig. ^, corre- 
sponding decay rates, v. Ma oc t" are in Table [| Note that 
the display is logarithmic, with an average Mach number 
exceeding 10 maintained for the first 100 yr in the low-speed 
case. Here we define the average Mach number as 



Ma - 



(cl> 



(11) 



where 



2 = (3.5-3/)(2.2-2/) e 
(3.3 - 2/) 2 p 

is an adiabatic speed of sound, and the bracket notation 
(...) means averaging over the domain: L~ 3 J2i j k- The 
Mach number decaying behaviour can be qualitatively ex- 
plained by noticing that it should be roughly proportional 
to the fraction {Ekin/ Eth) 1 ^ 2 , although the equality v — 
(77 — (J,) /2 is not generally true. 

The increase of the Mach number during the period 
« [30, 120] yr in the runs with rms velocity of 60 km s" 1 is 
due to the strong dissociative cooling at that time, as shown 
by the thermal energy decay curves, i.e. where \fj,\ > \n\. Note 
that this case, with the highest initial Mach number, pos- 
sesses the lowest Mach number once the shocks have formed. 
This is due to the dissociation, generating a warm atomic 
gas. The question then is: why do the turbulent motions 
decay so fast in this case, despite the low Mach number? 
The suggestion is that the shocks are preferentially running 
into the denser molecular gas, which can more efficiently 
dissipate the energy because it is cooler. In isothermal gas, 
however, this effect would not be apparent. 



3.2 Molecular hydrogen evolution 

The molecular fraction '/ ' is displayed in Fig. The three 
initial states correspond to three distinct physical regimes. 
With a rms velocity of 15 km s~ , very few dissociative 
shocks develop but some localised dissociation still occurs. 
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Figure 5. Histograms of the molecular fraction distribution. The black curve corresponds to the left y-axis and represents the volume 
distribution of molecules: the number of zones with value of molecular fraction within the interval [/, / + A/], A/ = 10 — 2 . The grey 
curve corresponds to the right y-axis and represents the mass distribution of molecules. The data were taken from the 128 3 run with 
rms velocity of 60 km s -1 . 



50 yr 



0.5 
0.4 
0.3 
0.2 



4.6e-23 8.0e-17 1.6e-1B 
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1 10 yr 
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density [g cm" 3 ] 



Figure 6. Histograms of the average molecular fraction distri- 
bution with density at different times. Molecules reform quickly 
even in regions with low densities. Data from 128 3 resolution run 
with initial rms velocity of 60 km s" 1 . 



With 30 km s -1 , a few per cent of the molecules are dis- 
sociated whereas at 60 km s _1 , the gas becomes over 80% 
atomic. 

Note that the minimum value of / is reached after 
~ 4ti. The dynamics alone determines the period of the 
dissociation phase, which does not overlap with the refor- 
mation phase for the chosen parameters. 

Reformation of molecular hydrogen is unexpectedly 
rapid. The expected H2 reformation time at 20 K and 



larger than the simulation time. At 100 yr, the temperature 
is ~ 80 K, predicting a reformation time of t R = (k R n ) _1 = 
2,000 yr. Yet, reformation is occurring over ~ 400 yr. This 
speed up is caused by the turbulence itself: the molecules 
preferentially reform in the denser and cooler locations. As 
weak shocks propagate through the gas, different regions 
are compressed and expanded. Hence the reformation time 
is not only controlled by the 'average' reformation time, but 
also by the strength of the turbulence. Given a turbulent 
dynamical timescale shorter than the average reformation 
timescale, then we can expect reformation to be accelerated. 

We have checked the degree of convergence of both the 
dissociation and reformation processes. Molecule reforma- 
tion is a gradual process and is, therefore, accurately rep- 
resented in the simulations. High resolution is, however, 
paramount to correctly follow the degree of dissociation. 
This is critical to the 30kms _1 turbulence since there are 
numerous intermediate speed shocks, partially dissociating 
the gas. For the low and high speed examples, however, the 
dissociation is basically zero and complete, respectively. Ex- 
haustive resolution studies of one-dimensional shocks with 



this code are presented elsewhere (Smith & Rosen 2002) 



10 b cm _-:1 is t R = {k R n )~ 



3,200 yr (Eq. 1) a factor of 5 



The remarkable manner in which the molecular fraction 
changes is illustrated in Fig. [5] for the high speed case. This 
figure displays the number of zones N z (/) with a molecular 
fraction in each interval [/, /+ A/] where A/ is taken as 
10" 2 . After 10 yr, there is an even distribution, suggesting 
that dissociative shocks have influenced about half the gas. 
Up to 50 years, the majority of zones have low values of /. 
Thereafter, the reformation has the following properties. 

• A single maximum in N z (f) shifts with time to higher 
values of /. 

• The maximum remains at N z ~ 10 5 zones. 

• The width of the peak is roughly constant (full width 
at 80% maximum is 0.16 ±0.02) 
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Figure 9. The correlation length changes with time in runs with 
different initial rms velocity. The data taken from the 128 3 run. 
Thick lines represent correlation length, thin line represent inverse 
average Mach number change. 



1 5 km s 




The autocorrelation function is denned as follows: 

L 

£(C(M)-<C>)(<(i + *,i + lO-(0) 
Sz(x,y) = — , 



(12) 



£(c(*,i)-<c» s 



where C,{x,y) is the quantity (integrated along the LOS): 



1 

CO,2/) = T^fel/^). 



(13) 



where £(x, y, z) is the density (or molecular density), and 
(£) is an average: 



(0 



L- 



(14) 



As a qualitative measure, however, it is more convenient to 
use an averaged correlation function, i.e. polar angle inte- 
grated: 



S 2 (r) 



tt/2 



S2 [r cos </>, r sin 4>)t d0, rE[0,L]. (15) 



Having integrated the correlation function, we have lost all 
information about anisotropic properties of the molecular 
density map. However, in our studies we have found that 
considerable deviations from polar symmetry of a correlation 
function occur usually only in regions of weak correlation, 
i.e. separated by a distance which is large in comparison with 
the correlation length. This suggests that large scale struc- 
ture which might seem to appear in molecular density maps 
(e.g. see the map on Fig. pi) is not coherent, and therefore 
turbulence in the simulations can be considered as isotropic. 

We find that the correlation length for the molecular 
column density grows as the turbulence decays (see Fig. ^ 
also compare Fig. | ^ and Fig. B ) as has been predicted for gas 



• The distribution of f hv mass disnlavs a single wide in general by |Mac Low (199$ and |Mac Low fc Osscnkopl 



peak ( grey line in Fig. 



(2000) ) . Defining the correlation length as the length over 



• Even in regions with a very low density, reformed 
molecules are found after 100 yr (Fig. |J). 

These vital results imply that / is dominated by a small 
range in values at any instant, at times after 100 yr. For 
this reason, three dimensional structure and maps based on 
column density of either all the gas or just the molecules 
appear almost identical. In this sense, isothermal simulations 
can be utilised to determine the structure of molecular clouds 
or cores but not the underlying fraction of molecules. 



which the average autocorrelation function has decayed by 
a factor of 10 (&(^ C r) = 0.1), we find that in the run with 
resolution of 128 3 , the correlation length is related to the 
average sound speed. The growth of £ cr shows a good corre- 
lation with the inverse average Mach number behaviour (see 
Fig. |^). An explanation for this is easily seen from the fol- 
lowing analysis. Defining the crossing time as t cross = L/v, 
where v is the average speed in the box, we find that the 
'information' during this time has travelled as much as 
£cr = c~s x t C ross- From this, one immediately obtains the 
following relation: 



3.3 Statistical analysis 

In recent years, several studies of probability density func- 
tions in numerical simulations have been advanced as a 



step in their full statistic a l characterisation (e.g. Vazquez 



Semadeni fc Garcia (2001 ); Burkert fe Mac Low (2001)). The 



goal is to quantitatively compare observed maps of the inter- 
stellar medium with the predictions of simulations. Here, we 
present the first analysis of the results for the column den- 
sity of molecules from turbulent simulations which include 
the molecular chemistry. 



If = * « Ma 
L v 



(16) 



The autocorrelation functions for density and molecular 
density maps indeed look very similar, as expected given 
the above results for the evolution of the molecular fraction 
within a relatively narrow range of values. 

An analysis of the probability density functions (PDF) 
reveals no particular difference between PDFs for density 
and molecular density maps and PDFs from different mo- 
ments of time. Linear fit coefficients n of the histograms of 
( in log-linear coordinates (N oc (C/(max) K ) are distributed 
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Figure 7. Autocorrelation functions for molecular density and snapshot of molecular density map (from 256 3 data) at 100 yr. The 
molecular density map is the result of the integration of the data cube (fraction X density) along z axis. 



averaged CF, 600 yr 



molecular density map 
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Figure 8. Autocorrelation functions for the molecular density and snapshot of the molecular density map (from 256 3 data) at 600 yr. 
The molecular density map is the result of the integration of the data cube (fraction X density) along z axis. 



in a quite wide range: k ps [4.0, 6.8 ] w ith a mean value of 
5.13 and dispersion « 0.66 (see Fig.|lO|). 

It has been proposed that observational me asurements 
are oft en confine d to within a correlation length (Burkert & 



Mac Lew (2001)). Linear fit coefficients of the average his- 



tograms (mean from a hundred histograms) of regions which 
have size of about a correlation length, yield a distribution 
with a smaller dispersion (m 0.22) and have a mean value of 
2.90 (Fig. [j^). This is roughly consistent with the observed 
value for molecular clouds found by Williams et al. (2000); 



Blitz & Williams (1997). In their studies they have used 



two-dimensional column density maps of clouds observed in 
different radial velocity bins for an optically thin molecular 
species. Each pixel in the data cube (galactic latitude, longi- 
tude, radial velocity) has an associated antenna temperature 
T a , proportional to the column density of gas in that pixel. 
The PDFs was then defined as the total number N of pix- 
els with a certain column density C,/C, m ax, (where C, ma x is 



the maximum column density found in the data cube). For 
molecular clouds in Taurus the following asymptotic^] value 
of k is reached: k — 2.7 ± 0.4. Although our simulations 
are not directly comparable with the observations of Giant 
Molecular Clouds, the consistency is quite remarkable. 



3.4 Spatial Distribution 

In Figs. and [l§ we present representations of the three di- 
mensional distributions. These figures demonstrate how the 
molecular density and total density distributions converge 
between 50 and 100 yr. Initially, of course, while thermal 
dissociation is occurring in the shock waves, there is less 
correlation. 



1 k is found to be sensitive to the spatial smoothing, for smooth- 
ing beam > 20' the PDF reaches an asymptotic value 
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Figure 10. Probability distribution functions. The left figure presents histograms from the 60 km s _1 run at lOOyr. The right figure 
presents histograms from the 30 km s _1 run at 100 yr. The thicker line corresponds to the averaged PDF from areas with a size of about 
or less then one correlation length (sa 0.23L), with the associated linear fit coefficient K cr . The data are taken from the 256 3 run. 



The general structure of the clouds is hard to define. 
There are some clumps and filament-like structures but little 
visual evidence for sheets. Some clumps and filaments are 
surrounded by diffuse envelopes, which are adjacent to other 
filaments and clumps. Self-gravity or sweeping by large-scale 
shocks could generate more coherent structures. 



4 CONCLUSIONS 
4.1 Summary 

We have presented the properties of a specific model for 
molecular turbulence. We carried out three dimensional hy- 
drodynamical simulations of decaying supersonic turbulence 
in molecular gas. We included a detailed cooling function, 
molecular hydrogen chemistry and equilibrium C and O 
chemistry. We studied three cases in which the applied ve- 
locity field straddles the value for which wholesale dissocia- 
tion of molecules occurs. The parameters chosen ensure that 
for the high-speed turbulence, the molecules are initially de- 
stroyed in shocks and gradually reform in a distinct phase. 
We find the following. 

• An extended phase of power-law kinetic energy decay, 
as in the isothermal case, after an initial phase of slow dis- 
sipation and shock formation. 

• The thermal energy, initially raised by the introduction 
of turbulence, decays only a little slower than the kinetic 
energy. 

• The reformation of hydrogen molecules, as the fast 
turbulence decays, is several times faster than expected 
from the average density. This is expected in a non-uniform 
medium, as a consequence of the volume reformation rate 
being proportional to the density squared. 

• The molecules reform into a pattern of filaments and 
small clumps, enveloped in diffuse structure. 



50 yr 




150 yr 



250 y 




500 y 



Figure 11. Molecular density distribution snapshots. Data taken 
from the 128 3 resolution run with rms velocity 60 km s _1 . At 
early stages (when average molecular fraction is low) the dis- 
tribution of molecular density is different from mass density (see 
Fig. h2|); later, when hydrogen reforms, the distributions look very 
similar 



• During reformation, the remaining turbulence redis- 
tributes the gas so that the fraction of molecules is dis- 
tributed relatively evenly. Hence, the density and molecular 
density are almost identically distributed at any one time 
after 100 yr. 
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50 yr 



150 yr 




250 yr 



Figure 12. Mass density distribution snapshots. Data taken from 
128 3 resolution run with rms velocity 60 km s — 1 . Compare the 
distributions with the molecular density distribution on Fig. |ll| 



• The correlation length in our simulations grows with 
time as has been predicted for a gas in general. 

• The probability density functions sampled from regions 
with size of about one correlation length are consistent with 
observations. 



4.2 Discussion 

We mainly wish here to emphasise the insight these simula- 
tions provide into how molecular chemistry and supersonic 
dynamics combine. We have found that isothermal simula- 
tions are indeed very useful, not only for the rate of energy 
decay but also to trace the molecules. 

A simple reason for the fast decay is that a sufficient 



numb er of strong shocks survive. As shown by Smith et al 



(2000), the rate of energy decay in decaying turbulence is 



dominated by the vast number of weak shocks. These shocks 
are less efficient at energy dissipation. A second possible rea- 
son is that the curved shock structures create small scale 
vorticity, which leads to enhanced dissipation of kinetic en- 
ergy. It is not clear, however, why relatively more vorticity 
should be created when stronger shocks are present. Simu- 
lations of isolated curved shocks would help clarify the dis- 
sipation paths. 

It is plausible that the high Mach number turbulence 
could create more small- scale structure, l eading to a faster 
decay rates as argued by Mac Low (1999 ). 

The simulations presented here may directly aid our 
interpretation of fast molecular shocks in dense regions, 
such as Herbig-Haro objects. For example, supersonic decay 
would app ear to occur in the w ake of bow shocks such as 
in DR 21 (Davis fc Smith 1996). These simulations suggest 



that the turbulence decays rapidly and bow shock wakes will 
be bright but short. 

The fast reformation of molecules indicates that 
molecules may form on shorter time scales than previously 
envisaged. Hence molecular clouds may appear out of atomic 
clouds several times faster than anticipated in non-turbulent 
scenarios. Recently, some evidence has been found for rapid 



cloud formation and dissipation ( Ballesteros-Paredes et al 



199£; Hartmann et al. 2001). This implies that much of the 



interstellar medium may be undergoing both rapid dynam- 
ical and chemical changes, driven by sources of supersonic 
turbulence. 

The simulations analysed here provide a basis for much 
further work. For example, we can now investigate the prop- 
erties of hydrodynamic clouds in which dust or chemical 
abundances are non-uniform, clouds in which the turbulence 
is driven uniformly or non-uniformly, and clouds which are 
initially atomic. 
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cloud evolution. These reactions are: 

H + Hi — >H 2 (Al) 

H 2 + H i — ► 3H (A2) 

H 2 + H 2 i — > 2H + H 2 (A3) 

O + H 2 i — > OH + H (A4) 

OH + Hi — >0 + H 2 (A5) 

OH + C i — > CO + C (A6) 

CO + H i — > OH + C (A7) 

OH + H 2 i — ► H 2 + H (A8) 

H 2 + Hi — ► OH + H 2 (A9) 

The first three reactions ( |Al| , |A2| , |A3| ) are included in a semi- 
implicit cooling-chemistry step to calculate the temperature 
and H 2 fraction. Reaction ( |Al| ) takes place on grain surfaces 
with the rate given by equation (see section 2.2). 

We fix free oxygen and carbon abundances, ??o(0) and 
T)o(G). We then calculate the equlibrium O, OH, and CO 
abundances, assuming that the CO reactions are much faster 
than the H 2 reactions. Then, the remaning free oxygen is 
distributed into O, OH, and H 2 0, according to equlibrium 
abundances. In this manner we can calculate quite accu- 
rately the chemistry within fast shock s without overloading 
the hydrocode. For further details (see Smith & Rosen 2002, 
Appendix B). 

The cooling fuction used in the numerical code (see 
equations |l] - ^ ) is composed of 11 separate parts (one 
of which heats the gas): 



A = E A * 



(A10) 



The components are summarised in Table Al below 



APPENDIX A: THE CHEMISITRY AND THE 
COOLING FUNCTION 

We consider only a limited network of chemical reac- 
tions, which has been tested in one-dimensional simulations 



(Smith fc Rosen 2002) and which is critical for molecular 
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Table Al. Components of the cooling function 
Formulae 



Description 



Ai = Ai X n , where, on assuming standard dust properties, 



Ai 



3.8 x 10" 



3 T°' 5 (T - T d 



t (1.0-0.8exp(-75/T)) [ergs" 1 ] 



gas-grain (dust) cooling (Hollenbach & McKee 1989). In 
present calculations we fix Td ust = 20K. 



A2 



r(h) 



rotational coefficients at Mgh and tow density are: 



where the vibrational and 



1.10 X 10 _is exp(-6744/T) [ergs" 



collisional cooling associated with vibrational and 
rotational modes of molecular hydrogen. These f or- 
mulae based on equations 7 - 12 in |^epp fc Shull (1983 ) 



4° =8.18 x 10- 13 exp(6840/T) (n H fc H (o,l) + "h^h^oj)) [ergs" 1 ], 

the terms fcu(C> ,1) an d ^H 2 (0,l) are u : — > 1 collisional excitation 
rates, 

_ J 1.4 x 10~ 13 exp(T/125 - T 2 /577 2 ), 

H(0,1) ~ {1.0 x io- 12 r°- 5 cx P (-iooo/r), 

where T v = 1635 K, 



if T < T v , 
if T > T„, 



«H 2 (0,1) 



1.45 x 10" 



'' exp(-28728/(T + 1190)); 



(h) 



dcx{-19.24 + 0.474a; - 1.247a; 2 }, if T < T r 



3.9 x 10- 19 exp(-6118/T), if T > T r , 

where T r = 1087 K, and x = log(T/10000), 

_ Jdex{-22.90 - 0.553 - 1.148a; 2 }, if T <T h 
Q(n) ~ 1 1.38 X 10' 22 exp(-9243/T), if T > T h 

where T; = 4031 K, and Q(n) = (nu 2 )°' 77 + 1.2 (n H ) ' 77 
A3 = (n H ) 2 x A 2 



A 4 = (n H2 + 1.39n H ) x n H2 o x A3; A3 = 1.32 x 10- 23 (T/1000) a , 
where a = 1.35 - log(T/1000). 



collisional cooling- of ato ms. We have used Table 10 of 
Sutherland fc Dopita (1993) (with Fc = -0.5) for the form 



of A2 and we added an extra thermal bremsstrahlung term 
equal to 1.42 x 10~ 27 T ' 5 for T > 10 4 K 
cooling through rotational modes of water induced 
by collisions with both atomic and molecular hydro- 
gen. Given values of a fits the values tabulated by 



A 5 = 1.03 x 10- 26 n H2 n H2 oTexp(-2325/T)exp(-47.5/T 1 /3) 



fc Kaufman (1993) 



Neufeld 



cooling through vibrational modes of w ater induced 
hv collision s with molecular hydrogen (Hollenbach 



A 6 = 7.40 x 10- 27 n H A 3 nH 2 orexp(-2325/T)cxp(-34.5/T 1 / 3 ) 



McKee 1989) 



cooling through vibrational modes of water induced 
by collis ions with atomic hydrogen ( Hollenbach & Mc 



A 7 = 7.18 x 10 -1 " {(n H2 yk D , H 2 +"H"H 2 fcr>,H) 1 where dissociation 
coefficients are taken to be, 

fc D , H = 1.2 X 10~ 9 exp(-52400/T) (0.0933exp(-17950/T)) /3 [cm 3 s -1 ]. 
k D ,H 2 = 1.3 X 10~ 9 exp(-53300/T) (0.0908 exp(-16200/T))' 3 [em's" 1 

13= [l.0 + n (2/(71^ - n" 1 ) -fn" 1 )]^ 

critical densities are fit by the following approximation, n\ = 
dex{4.0 - 0.416a; - 0.327a: 2 } [cm" 3 ] and ra 2 = dex{4.845 - 1.3a; + 
1.62a; 2 } [cm" 3 ]. 

A 8 = -A 7 £, where A 7 = k R (see Eq.0) and § = nra H (l-/3)7.18x 10~ 12 



Ag = ncQnkTavj'/ (l + n a /n cr + 1.5(n a /n cr ) ' 5 ) , where the mean 
speed of the molecules vp = ^/8fcT/(-7rmH 2 ) and 
10 6 (T/1000)°' 75 [cm" 2 ], a = 3.0 X 10- 16 (T/1000)-°' 25 [cm" 2 ]. Av- 
erage density, n a = 0.5(rtn + nn 2 \/2). 

A10 = 1.83 x 10- 26 riH 2 nco' T1 exp(-3080/T)exp(-68/T 1/3 ) 



Kee 1989) 



cooling from the dissociation of molecular hydrogen. 

Factor 7.18 X 10~ 12 [erg] is the 4.18 eV dissociation energy; 
ni - is the density critical for dissociation by collision of 
molecular hydrogen with atomic hydrogen, ?i2 - with itself. 



heating resulting from the reformation of molecu- 
lar hydrogen. We employ j3 to parametrise the fraction of 
released energy which is thermalised rather than radiated, 
cooling through rotational modes of carbon monox- 
ide induced by collision with both atomic and 
molecular h ydrogen. We bas ed our equations on 



An = 1-28 x 10- 24 n H n C oT°- 5 exp(-3080/T)exp(-(2000/T) 3 - 43 ) 



Eqs. 5.2-5.3 in [McKee et al. (1982 ) 

cooling through vibrational modes of carbon 
monoxide i nduced by collisions w ith molecular hy- 
drogen (see [Neufeld fc Kaufman 1993; ) 

cooling through vibrational modes of carbon 
monoxide induced by collisions with atomic hydro- 
gen 



